Functional and structural reorganization in brain tumors: a machine learning approach using desynchronized functional oscillations

Neuroimaging studies have allowed for non-invasive mapping of brain networks in brain tumors. Although tumor core and edema are easily identifiable using standard MRI acquisitions, imaging studies often neglect signals, structures, and functions within their presence. Therefore, both functional and diffusion signals, as well as their relationship with global patterns of connectivity reorganization, are poorly understood. Here, we explore the functional activity and the structure of white matter fibers considering the contribution of the whole tumor in a surgical context. First, we find intertwined alterations in the frequency domain of local and spatially distributed resting-state functional signals, potentially arising within the tumor. Second, we propose a fiber tracking pipeline capable of using anatomical information while still reconstructing bundles in tumoral and peritumoral tissue. Finally, using machine learning and healthy anatomical information, we predict structural rearrangement after surgery given the preoperative brain network. The generative model also disentangles complex patterns of connectivity reorganization for different types of tumors. Overall, we show the importance of carefully designing studies including MR signals within damaged brain tissues, as they exhibit and relate to non-trivial patterns of both structural and functional (dis-)connections or activity.

Functional MRI samples neural signal every a specific number of seconds.The time between two volume acquisitions, known as the repetition time (TR), is usually kept constant throughout the study.As already stated in the main text, the TR was accidentally changed from 2100 ms to 2400 ms.The effects of acquiring fMRI with different TRs have been studied for a considerable amount of time (see for example 1,2 ).Briefly, longer TRs offer better signalto-noise ratios (SNRs) and are not recommended for fixed time paradigms (i.e., scans of fixed duration).Furthermore, the number of volumes should be first selected and then be kept constant throughout the experiment.Therefore, resting-state data offers more reliable SNRs mainly due to the need of rapid sampling during task conditions 3 .The dataset we studied contained resting-state scans with relatively long TRs.Most importantly, after the accidental TR shift, the time of acquisition was not kept fixed (see Methods).Hence, it resulted in slightly longer series with the same number of volumes.
Nonetheless, the discrepancy between subjects needed to be addressed accordingly.Recall that we based our technique on discrete Fourier decompositions where the properties rely on the sampling frequencies and the sampling rates 1 .Given a TR and a signal of length , the difference between two consecutive Fourier sampling frequencies is a ratio which should be preserved even in the case of having two different sampling rates or TRs.This simple heuristic yields a relationship between the number of points of the two signals: the ratio between  and the TR should be constant.Therefore, series with shorter TRs were extended Δ points by 0-padding them, with where TR 1 > TR 2 and 1 extra point is added after ensuring that the division outputs an integer number.After computing the transformation, the resulting frequency spectrum did not share the same length.Consequently, the spectrum from the shorter series was simply truncated at the same length as the rest of the spectrums.After inspecting the original spectrums (see Fig. 1; see also Fig. S2-S5) we noticed that no power was allocated to high frequencies thus ensuring no power loss after the truncation.We provide a comprehensive example with a sinusoidal signal consisting of only two frequencies in Fig. S14.Noteworthy, the total power of the signal after the extension is not kept constant, but the distribution of the power is correctly maintained (Fig. S14).Since the DAS defined in the text does not care about absolute power but rather their distributions, this approach safely preserved information of the series' dynamics.However, all the absolute power analyses conducted in this study were done considering the unpadded signals -although the statistical results were identical to those ones obtained with the padded and truncated signals (i.e., up to the third decimal point).

Generation of post-surgery graphs with artificial neural networks
A connectome for the -th subject at a particular time point  which, in our case, was mapped to control, pre-and post-surgery stages   ,   and   .Only the lower triangular part of the connectomes was flattened into a vector   () of  components, given that structural connections are always symmetric and no self-loops were considered.Each one of these components   () with  = 1, … ,  represents a link of strength  between two brain regions.We simplify the notation and refer to pre-surgery connectomes as   , to postsurgery connectomes as   and to healthy connectomes as   .
Following previous work 4 , we defined a distribution of binary links based on connectomes from   healthy subjects.For the -th edge,   = 0,1 is a Bernoulli distributed binary variable with probability where Θ(⋅) is the Heaviside step function and  = 0.2 is a threshold that filters links with insufficient strength and therefore minimizing the false positive rate.For the choice of the threshold, we considered the Shannon entropy and the modularity 5,6 of the network built according to Eq. (S2).We observed that small thresholds maintained a high entropy allowing for the inclusion of more connections without excessively altering the structure of the prior (Fig. S7).Lower thresholds allow more variability in the generated prior since more spurious connections can be considered.
For each patient and edge, the post-surgery probability of having a meaningful link   = 1 with strength  is conditioned on the pre-operative graph.We exploit this fact with the wellknown Bayes' theorem, is the likelihood of having a post-surgery -strengthened connection conditioned on the pre-surgery connectome.For Eq. (S3) to hold, (λ  = 1) must be independent on the pre-surgery graph   .This assumption is reasonable since the anatomical prior was built using only healthy controls, therefore, not considering lesioned connections.We sampled each -th connection using the Maximum A Posteriori criterion to generate the post-surgery graphs.To train the network, the sampled connections were then compared to the ground truth using the Mean Squared Error (MSE).
Although many possibilities emerge for estimating the likelihood in Eq. (S3), we used a fully connected network with one hidden layer of size /2.Therefore, for each edge  in the -th connectome, the likelihood in Eq. (S3) was obtained by propagating the input through the artificial network as follows: where () = 1 1+ − ∈ [0,1] is a sigmoid activation function.Note how the network compresses the information in the intermediate layer forcing a low dimensional representation of connectomes' features.The model weights   input and   output were updated using stochastic gradient descent to minimize the mean squared error between the connectomes generated by the model and the ones obtained by applying the tractography pipeline described in the main text.Several architectures were considered, but we found that adding more layers or even considering 1D or 2D convolutions did not add significant improvements to the model.We report the mean results for a two hidden layer fully connected layer in Table S2.

Weight probability distribution of brain networks
For all graphs, we computed the degree of each node by summing over all the connections where  is the number of nodes, 166 in our graphs.To compute the probability distribution of a node having degree , we counted all the nodes whose degree fell in a bin of width Δω as where  is again the number of nodes that acts as a normalization constant.Networks wired at random have a binomial degree distribution, however it was shown that biological networks deviate from this tendency towards a scale-free or, in the specific case of brain networks, a broad-scale distributions 7,8 .In both cases, the number of highly connected nodes (i.e., nodes with large degree) decays following a power-law in the case of scale-free and an exponentially truncated power-law for broad-scale distributions.However, for weighted networks, it has been shown that a lognormal distribution is better suited to describe brain connections 9,10 .The probability density of a lognormally distributed variable follows where ,  are the geometric mean and standard deviation.Note, that the main manuscript deals with the variable 1 + , therefore Eq. (S6) should be corrected for this.Given that the Jacobian of the transformation is 1, the proof is straightforward.
Therefore, the variable 1 +  also follows a lognormal distribution as seen in Fig. 6 of the main text.We adjusted a power law distribution to the weighted degrees 11 , where  is the power law exponent measuring the decay in the probability of finding strongly connected nodes in the network (e.g., hubs).We independently fit a power law distribution to each network derived from the healthy cohort and patient.Importantly, for each patient, we fit a power law distribution for the network obtained with our pipeline and another fit for the network obtained with a canonical multi-shell multi-tissue approach without masking the tumor.Next, we compared the differences between the exponent  for each pipeline and the healthy cohort and the differences between pipelines (see Fig. S6d).Importantly, we consider two different scaling ranges:   = 100 and   = 300, the latter one evaluating the scale-freeness in the large asymptotic range.S3 Leave One Out cross validation results and a multi-shell multi-tissue fiber reconstruction pipeline.Results for the tested subject in each trained model when using only a multi-shell multi-tissue fiber reconstruction method.(MSE: Mean Squared Error; MAE: Mean Absolut Error; PCC: Pearson Correlation Coefficient; CS: Cosine Similarity; KL Kullback-Leiber and JS: Jensen-Shannon Divergences).For each fold the subject tested was not included in the training set (see Methods).The best subject in each according to each metric is highlighted in bold.Correlations between all pairs of metrics were high (|r|>0.8and p<0.01).Only in PCC and CS a higher number means a better fit.Bold numbers indicate better performance; lower is better for MSE, MAE, KL, and JS while higher is better for CS and PCC.

Quantifying the effects of the Repetition Time (TR) on the reorganization of the functional networks
To further understand the effects that the inconsistent TRs might have introduced to the analyses carried in this work we studied the same scores after excluding the 7 patients and 4 controls with TR=2.1s.While these types of procedures may effectively reduce the statistical power of the analyses, they can help to elucidate hidden confounders in the data.The reduced group showed a medium correlation between the network similarity measure and the DAS of the DMN (r=-0.215,p=0.392, df=18, two-tailed exact test) and a strong correlation between the alterations in the lesion and DMN (r=-0.603,p=0.008, df=18, twotailed exact test).
We also did permutation analyses to determine whether the 7 excluded subjects with brain tumors represented a unique subsample or rather a random instantiation of 7 patients.We computed  = 1000 correlation coefficients between the same three magnitudes mentioned in the previous paragraph and in the main text after randomly excluding 7 patients from the pool.The exclusion was done by resampling without replacement from the patients' ID.For the network similarity analyses the effect of the TR appeared to be significant although it didn't alter the interpretation nor the conclusion of the main analyses (Fig. S15a).For the relationship between local and global alterations in the DMN, the TR did not have a significant effect (Fig. S15b).
To exclude the possibility of reporting false inferences, we measured the exact contribution of the TR in network similarity measure.For that, we fitted two linear models that explained the differences in the Pearson correlation between the networks of the patients and controls as a function of the DAS and the TR: Node similarity =  ⋅  + intercept, (9) Node similarity =  ⋅  +  ⋅  + intercept.
(10) Importantly, this analysis was performed with the whole dataset.For the first model in Eq.S9, which only included the DAS, the effect was, as expected, significant (a=-0.057,p=0.01, df=23, two-tailed t-test).The coefficient of determination was identical to the one in Fig. 2a; again, as expected.However, when explicitly including the TR, the coefficient in Eq.S10 quantifying the contribution of the TR to the node similarity measure was non-significant (b=0.178,p=0.139, df=22, two-tailed exact-test, CI=[-0.062,0.418]).
These results can be interpreted in the following way.The padding strategy correctly removed all the effects of the TR for the analyses that were performed solely on the frequency domain (e.g., local and global alterations).Yet, for the other case, although the contribution of the TR to the results was proven to be non-critical, the TR inconsistency remained slightly observable.However, and most importantly, all the interpretations and conclusions derived from the analyses in the main text remained identical and only decreased in size, as was expected given the reduced statistical power (Fig. S15c-f).
11. Alstott, J., Bullmore, E. & Plenz, D., powerlaw: A Python Package for Analysis of Heavy-Tailed Distributions.PLoS ONE 9, e85777 (2014).for the control population, while "*" codes for statistical significance.The relative power was obtained by averaging over subjects and DMN regions and taking healthy subjects as the baseline (i.e., equal to 1).b Functional DMN from the same patient and mean of control subjects (see Methods).c Functional complexity as measured by the distribution of correlations for the patient's network and the mean of control [mean ± SEM] (see Methods).d Functional signal of two randomly selected regions from the same patient (red) as opposed to all the control subjects (light blue).No apparent difference is found between raw time series across regions, subjects, and patients.Even if the tumor is not masked, the algorithm damps the anisotropic component of the diffusion signal.c Tumoral and peritumoral FODs reconstructed with the SS3T algorithm applied only inside the lesion mask.Only this algorithm can recover the diffusion signal, which in this example, loosely maps the inflammation caused by the tumor.d Differences in the scale-freeness, as measured by the power law exponent, between the networks obtained with our pipeline and a canonical multi-shell multi-tissue approach (blue scatter plot) and between each pipeline and the healthy cohort (black boxplot).Green triangles depict the means, black horizontal lines depict the medians, and the boxes show the 1 st and 3 rd quartiles.
Significance tests in green show the p-values for the differences between means.  refers to the minimum scaling range taken to fit the power law distribution.The higher the value, the more extremal values are considered.Only the absolute values were significantly different from zero (oneside t-test and U-test).d Scatter plot of the two components found via Fast Independent Component Analysis applied to the DAS between patients and control subjects.The same two clusters (red, and blue) were consistently found.e TOP, Alteration in dynamics between different groups of tumors.BOTTOM, Logistic regression between DAS and periventricular (PV) tumor patients.Qualitative inspection reveals a tendency for periventricular tumors to display slower BOLD dynamics, although they were found to be non-significant presumably due to the small number of samples (N=5, p=0.33, one-sided U-test).Only the absolute values were significantly different from zero (one-side t-test and U-test).d Scatter plot of the two components found via Fast Independent Component Analysis applied to the DAS between patients and control subjects.The same two clusters (red, and blue) were consistently found.e TOP, Alteration in dynamics between different groups of tumors.BOTTOM, Logistic regression between DAS and periventricular (PV) tumor patients.Qualitative inspection reveals a tendency for periventricular tumors to display slower BOLD dynamics, although they were found to be non-significant presumably due to the small number of samples (N=5, p=0.28, one-sided U-test).The discrete Fourier sampling frequencies for the original series in A and the series after adding zeropadding.Note how the black and red hyphened lines coincide after the truncation (frequencies to the right of the gray vertical line are discarded).c, d Power spectrum and cumulative power of the two series before (c) and after (d) zero-padding.Note that the cumulative power in the short unpadded series reaches 100% before the long series.This results in artificially altered DAS.After padding and truncation, these effects are successfully corrected.).e Scatter plot showing the null correlation of DAS with distance and overlap between the lesion and the DMN centroids (only for subjects with TR=2.4s).f A strong linear trend was found between alterations inside the patients' tumor and alterations in the DMN as measured by DAS (only for subjects with TR=2.4s).The orange line corresponds to the significant linear fit (two-sided exact test), and the shaded areas mark the 95% confidence interval.

Fig. S1 :
Fig. S1: DAS example of 3 sinusoidal signals.a Two sinusoidal signals with the same frequency and in complete dephasing.Their temporal correlation is exactly -1 and their mean is always zero and no coupled behavior can emerge.b The power spectrum, however, is identical (completely overlapping).Both signals, therefore, carry the same amount of power in the same frequency band.The lines represent the cumulative power distribution as defined in Eq. (3).In short, these two signals are indistinguishable since the dephasing might be a consequence of different starting time points.c Two sinusoidal signals with different frequencies of oscillations and arbitrary dephasing can carry the same power.The mean of both signals is now the composition of both frequencies.d The spectrum of frequencies is, nonetheless, very different.Conventional distance metrics would only capture differences in the information content, but the Dynamics Alteration Score (DAS) can distinguish which one of the signals displays faster or slower oscillations.The values of both the power and the DAS are entirely conditioned to the scale used (i.e.,  ∈ [0,1] if we had not defined the distributions as percentages).

Fig. S2 :
Fig. S2: sub-PAT15 DMN. a Cumulative power (lines), power distribution (bars), autocorrelation function (ACF), and total power relative to healthy subjects (histogram) for the patient (red) and mean of controls (blue) functional signals of the regions in the DMN.Error bars [mean ± SEM] indicate the results for the control population, while "*" codes for statistical significance.The relative power was obtained by averaging over subjects and DMN regions and taking healthy subjects as the baseline (i.e., equal to 1).b Functional DMN from the same patient and mean of control subjects (see Methods).c Functional complexity as measured by the distribution of correlations for the patient's network and the mean of control [mean ± SEM] (see Methods).d Functional signal of two randomly selected regions from the same patient (red) as opposed to all the control subjects (light blue).No apparent difference is found between raw time series across regions, subjects, and patients.

Fig. S3 :
Fig.S3: sub-PAT24 DMN. a Cumulative power (lines), power distribution (bars), autocorrelation function (ACF), and total power relative to healthy subjects (histogram) for the patient (red) and mean of controls (blue) functional signals of the regions in the DMN.Error bars [mean ± SEM] indicate the results for the control population, while "*" codes for statistical significance.The relative power was obtained by averaging over subjects and DMN regions and taking healthy subjects as the baseline (i.e., equal to 1).b Functional DMN from the same patient and mean of control subjects (see Methods).c Functional complexity as measured by the distribution of correlations for the patient's network and the mean of control [mean ± SEM] (see Methods).d Functional signal of two randomly selected regions from the same patient (red) as opposed to all the control subjects (light blue).No apparent difference is found between raw time series across regions, subjects, and patients.

Fig. S4
Fig. S4 sub-PAT19 DMN. a Cumulative power (lines), power distribution (bars), autocorrelation function (ACF), and total power relative to healthy subjects (histogram) for the patient (red) and mean of controls (blue) functional signals of the regions in the DMN.Error bars [mean ± SEM] indicate the results for the control population, while "*" codes for statistical significance.The relative power was obtained by averaging over subjects and DMN regions and taking healthy subjects as the baseline (i.e., equal to 1).b Functional DMN from the same patient and mean of control subjects (see Methods).c Functional complexity as measured by the distribution of correlations for the patient's network and the mean of control [mean ± SEM] (see Methods).d Functional signal of two randomly selected regions from the same patient (red) as opposed to all the control subjects (light blue).No apparent difference is found between raw time series across regions, subjects, and patients.

Fig. S5 :
Fig. S5: sub-PAT02 DMN. a Cumulative power (lines), power distribution (bars), autocorrelation function (ACF), and total power relative to healthy subjects (histogram) for the patient (red) and mean of controls (blue) functional signals of the regions in the DMN.Error bars [mean ± SEM] indicate the resultsfor the control population, while "*" codes for statistical significance.The relative power was obtained by averaging over subjects and DMN regions and taking healthy subjects as the baseline (i.e., equal to 1).b Functional DMN from the same patient and mean of control subjects (see Methods).c Functional complexity as measured by the distribution of correlations for the patient's network and the mean of control [mean ± SEM] (see Methods).d Functional signal of two randomly selected regions from the same patient (red) as opposed to all the control subjects (light blue).No apparent difference is found between raw time series across regions, subjects, and patients.

Fig. S6 :
Fig.S6: Peritumoral anatomy and FODs.a Binary mask depicting the gray-matter (LEFT), whitematter (CENTER), and tumor (RIGHT) tissues of sub-PAT16.The tissue surrounding the lesion is automatically labeled as gray-matter thus imposing an anatomical barrier to regular fiber tracking algorithms.b Tumoral and peritumoral FODs reconstructed with the MSMT algorithm when explicitly masking the tumor (LEFT) and without masking the tumor (RIGHT).Even if the tumor is not masked, the algorithm damps the anisotropic component of the diffusion signal.c Tumoral and peritumoral FODs reconstructed with the SS3T algorithm applied only inside the lesion mask.Only this algorithm can recover the diffusion signal, which in this example, loosely maps the inflammation caused by the tumor.d Differences in the scale-freeness, as measured by the power law exponent, between the networks obtained with our pipeline and a canonical multi-shell multi-tissue approach (blue scatter plot) and between each pipeline and the healthy cohort (black boxplot).Green triangles depict the means, black horizontal lines depict the medians, and the boxes show the 1 st and 3 rd quartiles.

Fig. S7 :
Fig.S7: Choosing the threshold for the Anatomical Prior.a Entropy and Modularity as a function of the selected threshold.These two properties were rather robust in certain intervals of the threshold value.Therefore, the choice of the threshold inside each one of these intervals was not critical to the performance of FCNET.Vertical dashed lines mark approximately the frontiers of the aforementioned intervals.b Number of occurrences of the variable described by Eq.S2 for each one of the thresholds marked previously.Higher thresholds deleted a considerably large number of possible connections limiting the generalization power of FCNET.

Fig. S8 :
Fig. S8: DMN and tumor core BOLD signals.a Scatter plot showing the null correlation of DAS with distance and overlap between the tumor core and the DMN centroids.b A strong linear trend was found between alterations inside the patients' tumor core and alterations in the DMN as measured by the DAS.The orange line corresponds to the significant linear fit (two-sided exact test), and shaded areas mark the 95% confidence interval.c Differences in total power relative to healthy subjects (LEFT) and dynamics (RIGHT).Only the absolute values were significantly different from zero (oneside t-test and U-test).d Scatter plot of the two components found via Fast Independent Component Analysis applied to the DAS between patients and control subjects.The same two clusters (red, and blue) were consistently found.e TOP, Alteration in dynamics between different groups of tumors.BOTTOM, Logistic regression between DAS and periventricular (PV) tumor patients.Qualitative inspection reveals a tendency for periventricular tumors to display slower BOLD dynamics, although they were found to be non-significant presumably due to the small number of samples (N=5, p=0.33, one-sided U-test).

Fig. S9 :
Fig. S9: DMN and non-necrotic tumor BOLD signals.a Scatter plot showing the null correlation of DAS with distance and overlap between the non-necrotic tissue and the DMN centroids.b A strong linear trend was found between alterations inside the patients' non-necrotic tissue and alterations in the DMN as measured by the DAS.The orange line corresponds to the significant linear fit (two-sided exact test), and shaded areas mark the 95% confidence interval.c Differences in total power relative to healthy subjects (LEFT) and dynamics (RIGHT).Only the absolute values were significantly different from zero (one-side t-test and U-test).d Scatter plot of the two components found via Fast Independent Component Analysis applied to the DAS between patients and control subjects.The same two clusters (red, and blue) were consistently found.e TOP, Alteration in dynamics between different groups of tumors.BOTTOM, Logistic regression between DAS and periventricular (PV) tumor patients.Qualitative inspection reveals a tendency for periventricular tumors to display slower BOLD dynamics, although they were found to be non-significant presumably due to the small number of samples (N=5, p=0.28, one-sided U-test).

Fig. S14
Fig.S14Solving the TR inconsistency using Fourier analysis.a Two series sampled with different rates (TR=2400ms blue -TR=2100ms orange) while keeping the number of points constant at 180 points.As a result, the short TR series is shorter in time.The underlying signal consists of a combination of two sinusoidal signals, each one with different frequencies.Although the signal is the same, the sampled series differ slightly.Importantly, the TR mismatch in the real data is the same.b The discrete Fourier sampling frequencies for the original series in A and the series after adding zeropadding.Note how the black and red hyphened lines coincide after the truncation (frequencies to the right of the gray vertical line are discarded).c, d Power spectrum and cumulative power of the two series before (c) and after (d) zero-padding.Note that the cumulative power in the short unpadded series reaches 100% before the long series.This results in artificially altered DAS.After padding and truncation, these effects are successfully corrected.

Fig. S15 :
Fig. S15: Contribution of the Repetition Time to the reorganization of the network.a Null distribution of the correlation coefficients between the PCC and DAS of the DMN.b Null distribution of the correlation coefficients between the DAS of the DMN and the lesion.c Correlations between node similarity measured by the Pearson correlation coefficient and the Dynamics Alteration Score (DAS) of the Default Mode Network (DMN; only for subjects with TR=2.4s).The inset shows how alterations in the dynamics shape the organization of the network regardless of the direction of the displacement of the frequency spectrum.d Alterations in dynamics correlate with alterations in network complexity.Change in the complexity of the patients' DMN with respect to the healthy pool (Directionality ΔΘ and Magnitude |ΔΘ|) as a function of the DAS (only for subjects with TR=2.4s).e Scatter plot showing the null correlation of DAS with distance and overlap between the lesion and the DMN centroids (only for subjects with TR=2.4s).f A strong linear trend was found between alterations inside the patients' tumor and alterations in the DMN as measured by DAS (only for subjects with TR=2.4s).The orange line corresponds to the significant linear fit (two-sided exact test), and the shaded areas mark the 95% confidence interval.

Table S1 Leave One Out cross validation results and a hybrid fiber reconstruction pipeline.
Results for the tested subject in each trained model.(MSE: Mean Squared Error; MAE: Mean Absolut Error; PCC: Pearson Correlation Coefficient; CS: Cosine Similarity; KL Kullback-Leiber and JS: Jensen-Shannon Divergences).For each fold the subject tested was not included in the training set (see Methods).The best subject in each according to each metric is highlighted in bold.Correlations between all pairs of metrics were high (|r|>0.8and p<0.01).Only in PCC and CS a higher number means a better fit.Bold numbers indicate better performance; lower is better for MSE, MAE, KL, and JS while higher is better for CS and PCC.

Table S2 Comparison of alternative models (mean
Bold numbers indicate better performance; lower is better for MSE, MAE, KL, and JS while higher is better for CS and PCC. ± SEM).Comparisson of alternative models using a Leave One Out cross validation scheme in 6 metrics (MSE: Mean Squared Error; MAE: Mean Absolut Error; PCC: Pearson Correlation Coefficient; CS: Cosine Similarity; KL Kullback-Leiber and JS: Jensen-Shannon Divergences).The last row corresponds to the same FCNET model but using as inputs brain networks obtained with the alternative multi-shell multi-tissue fiber reconstruction pipeline (see Main Text; see also Methods).